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Abstract — This paper discusses an algorithm for estimating 
the safe maneuvering envelope of damaged aircraft. The al- 
gorithm performs a robust reachability analysis through an 
optimal control formulation while making use of time scale 
separation and taking into account uncertainties in the aerody- 
namic derivatives. Starting with an optimal control formulation, 
the optimization problem can be rewritten as a Hamilton- 
Jacobi-Bellman equation. This equation can be solved by level 
set methods. This approach has been applied on an aircraft 
example involving structural airframe damage. Monte Carlo 
validation tests have confirmed that this approach is successful 
in estimating the safe maneuvering envelope for damaged 
aircraft. 

I. INTRODUCTION 

All transportation systems need to focus on safety, but 
this applies especially to civil aviation. Therefore, in civil 
aviation, many developments focus on improving safety 
levels and reducing the risks that critical failures occur. In a 
recent study by CAST/ICAO, it can be observed that “loss 
of control in flight” (LOC-I) is the most frequent primary 
accident cause. This study is based on a statistical analysis 
of aircraft accidents between 2002 and 2011, and indicates 
that this category accounts for as much as 23% of all fatal 
aircraft accidents and involves most fatalities[l]. Benefit can 
be gained by developing technology which prevents these 
LOC-I accidents. From a flight dynamics point of view, 
with the technology and computing power available on this 
moment, it might have been possible to recover some of 
the aircraft in the accident category described above on 
the condition that non-conventional control strategies would 
have been applied. These non-conventional control strategies 
involve the so-called concept of fault tolerant flight control 
(FTFC), where the control system is capable of detecting and 
adapting to changes in the aircraft behaviour. 

One FTFC strategy option is using a model based control 
routine. Previous research focused on a physical model 
approach[2]. In this setup, experiments have shown that not 

*This work is supported by a Marie Curie International Outgoing Fel- 
lowship (IOF) within the 7th European Community Framework Program. 

1 T. Lombaerts is with the German Aerospace Center DLR, Institute of 
System Dynamics and Control, Department of Aircraft Systems Dynamics, 
82234 Wefiling - Oberpfaffenhofen, Germany. Currently, he is Marie Curie 
Fellow and visiting researcher at NASA Ames Research Center, Intelligent 
Systems Division, Adaptive Control and Evolvable Systems (ACES) Group, 
Moffett Field, CA94035, USA. thomas . lombaerts@dlr.de and 
thomas . lombaerts@nasa . gov 

2 S. Schuet is with NASA Ames Research Center, Intelligent Systems 
Division, Physics based Methods Group, Moffett Field, CA94035, USA. 

3 K. Wheeler and D. Acosta and J. Kaneshige are with NASA Ames Re- 
search Center, Intelligent Systems Division, Adaptive Control and Evolvable 
Systems (ACES) Group, Moffett Field, CA94035, USA. 


only a reconfiguring controller is needed, but also some form 
of flight envelope protection, which prevents the airplane 
from leaving the safe flight envelope and losing control in 
flight[2]. The main challenge in this context is determining 
the new bounds of the safe flight envelope after failure, which 
are then used by the envelope protection algorithm[3]. 

Alternative and complementary research approaches for 
the purpose of loss of control prevention and prediction 
are among others passive adaptive control[4], data-based 
predictive control[5] and real-time optimal envelope limit 
estimation[6]. 

Determination of the flight envelope has been done in 
the literature through various methods. The most straight- 
forward methods include wind tunnel testing, flight test 
experiments and high-fidelity model-based computation of 
attainable equilibrium sets or achievable trim points[7], [8], 
[9]. More complex methods include formulating flight enve- 
lope estimation as a reachability problem and solving this 
with level set methods and Hamilton- Jacobi equations [10], 
[11], [13], [14], [15], possibly with time scale separation [16] 
or semi-Lagrangian level sets [18]. Alternative methods rely 
on linearization and region of attraction analysis [19], deter- 
mining controllability/maneuverability limits in a quaternion- 
based control architecture [20] or robustness analysis for 
determination of reliable flight regimes [21]. An approach 
suggested by Boeing, as part of the NASA program Dy- 
namic Flight Envelope Assessment and Prediction (DFEAP), 
uses Control-Centric Modeling, dynamic flexible structure 
and load models [22]. In the frequency domain, stability 
margins can be estimated in real time via nonparametric 
system identification [23]. More focused techniques inspired 
by flight dynamics exist as well, such as determining the 
minimum lateral control speed[24]. 

From the perspective of the physical approach, the pre- 
ferred interpretation of the safe maneuvering envelope con- 
siders reachability from the trim envelope. The stable and 
controllable trim envelope is considered an a-priori safe set. 
The backwards reachable set is defined as the set of states 
from where the trim envelope can be reached. The forwards 
reachable set is the set of states which can be reached from 
the trim envelope. Then the safe maneuvering flight envelope 
is the cross section of the forwards and backwards reachable 
sets. This interpretation is illustrated in Fig. 1. In addition, 
the backwards reachable set is the survivable flight envelope. 
After an upset due to damage, turbulence, a wake encounter 
etc., one can bring the aircraft back to a safe trim condition 
if the current flight condition is situated inside the backwards 
reachable set. 



backwards 


operating set = safe flight envelope 

forwards reachable set 


Fig. 1. Safe maneuvering envelope as intersection between forwards and 
backwards reachability, modified from source: van Oort[18] 


The aim is to perform a combined forward and backward 
reachability analysis from the trim envelope as efficiently 
as possible, for on-line implementations. Based on previous 
research[6], level set methods are an excellent candidate. 
Two of the major challenges are the computational load and 
how to cope with nonlinear systems with higher dimensions. 
In general, an increase in technology readiness level (TRL) 
is envisaged. 

Nonlinear systems with higher dimensions can be simpli- 
fied by considering the principle of time scale separation[16]. 
The structure of time scale separation is analogous as applied 
for the fault tolerant control algorithm developed earlier[2]. 
The overview can be found in Fig. 2, which illustrates that 
a nine dimensional nonlinear problem is decoupled in three 
consecutive three dimensional optimization problems. 
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interval [t,t f ] to the set of admissible input values U. Define 
0(T,f,x,u(-),A) as the state trajectory. A are defined as 
parameter uncertainties inside the predefined set of expected 
uncertainty values D. Given a set of states K C W 1 , the 
reachability question can be naturally formulated regarding 
the relation between the set K and the state trajectories 0 
of Eq. (1) over the horizon T. Problem of interest is the 
following: 

Robust reachability: Does there exist auG^^j and a 
t G [0, T] such that the trajectory 0 of the state x satisfies 
x e K, irrespective of A? 

The optimization problem can be formulated as a pursuit 
evasion game over the horizon T > 0 with target set K CW 1 
[17]. It is assumed that u is trying to bring or keep the 
state in the set K , whereas A is trying to drive it out of 
K. To ensure the game is well-posed, u is restricted to 
play non-anticipative strategies with respect to the unknown 
uncertainties A. 

For the types of safety problems considered here, a set of 
initial states has to be established such that u can win the 
game, in other words the set & can be characterized: 

= { x G R" | VA e £>, 3u 6 &\tj] , 

3tg [t,T],<j>(T;,t,x,u(-),A)eK} 

As done elsewhere in the literature [10], the characteriza- 
tion of this set can be done according to the principle of 
duality: 

@(t,K) = (S(t,K c )) c (3) 

where # c stands for the complement of •. Through this 
principle, it can be characterized as an INFMIN problem[10]. 
The crux is to include the A’s as disturbances in the optimiza- 
tion function, they oppose the optimization over u. Consider 
a closed set K, that can be written as the level set of a 
continuous function / : W 1 — >> M, i.e. ^={xgR”|/(x) ^ 0}. 
As a consequence, the Invariance optimization formulation 
becomes[16]: 


Fig. 2. Separation of dynamics over high bandwidth, middle range and 
low bandwidth. F Ax and Fa z are defined in Eq. (9) and (10). 

with: 


J (t,K) = { x 6 R"| Vi (x,r) > 0} 


( 4 ) 


II. OPTIMAL CONTROL FORMULATION 

It has been shown in the literature that maneuvering 
envelope estimation through reachability can be reformulated 
in the optimal control framework[10]. Consider a continuous 
time control system: 


V\ (x,f) = inf sup min / (0 (r,fx,u (•) , A)) (5) 

Ae ■D T ^[ t : T ] 

This can be reformulated into an Hamilton-Jacobi-Bellmann 
Partial Differential Equation (HJB PDE)[10], [17]: 


x = f(x,u,A) (1) 

with x G R", u G U C R m , A G D C R*, f (•, •) : R" x U — > R", 

a function: 

I(-):R n ->R ( 2 ) 

and an arbitrary time horizon T > 0. Let W [^/] denote the 
set of Lebesgue and bounded measurable functions from the 


inf sup ^- l - (x, t) f (x, u, A) > =0 

u (‘) e ^,r] AeD J 

( 6 ) 

where V\ (x, T) = /(x) holds for backward integration and 
V\ (x, t) = l (x) applies to forward integration. These HJB 
PDE’s can be solved by level sets, for which a toolbox is 
available in Matlab® [11]. 


dV\ 

(x,C+ min 

dt K J t e[t,T] 


III. APPROACH 

The approach to calculate the safe maneuvering envelope 
after damage is based on the following steps: 

• Identify the updated aircraft parameters after damage 
(not discussed here, see [12]). This concerns primarily 
estimating new post-damage values for the aerodynamic 
derivatives such as C Lo , C La , C Do , C Da , C D ^ 2 and C 7/3 
together with their uncertainty bounds. 

• Calculate the post-damage trim envelope based on the 
updated aircraft parameters (not discussed here, see 
[ 12 ]). 

• Based on the previous step, define reference trim bound- 
aries for airspeed V, flight path angle y and bank angle 
(p as well as grid step size AV, Ay and A (p. 

• Define an implicit function accordingly over V and y. 
This needs to be done for every value of (p in case 
speed and flight path angle are bank angle dependent, 
i.e. V =V((p) and y= y(<p). 

• The Level sets method toolbox[ll] relies on the Hamil- 
tonian in Eq. (6) as a gradient to evolve the implicit 
function and thus reference boundaries over time 

• The cost function V\ becomes a three dimensional 
functions, where cross sections reflect the situation for 
a specific time instant t a . 

• A dissipation function is needed to guarantee numerical 
stability during these calculations. As a consequence, 
slightly conservative results are obtained for the bound- 
aries, but analysis has shown that this dissipation has 
a minor effect on the results. In this specific context, 
the chosen dissipation function is a Lax Friedrichs 
dissipation function[ll]. 

IV. APPLICATION EXAMPLE 

To illustrate how the envelope estimating algorithm works, 
a nonlinear 3D aircraft example is considered. At this point, 
only the slow dynamics as specified in Fig. 2 are considered. 
Future work will extend to the faster dynamics. The data used 
in this example are based on the RCAM (Research Civil 
Aircraft Model) simulation model[25]. The acting forces on 
the aircraft are illustrated in Fig. 3 for a symmetric flight 
condition. 





Fig. 3. Acting forces on the aircraft model, source: Lygeros[10] 


For the complete 3D situation, the equations of motion for 
V and y are written as follows [2]: 

Fj x — W siny — mV (7) 

7^ z cos(p+7^ F sin(p + lFcosy = —mVy (8) 

Where the aerodynamic force components can be simplified 
assuming small aerodynamic angles a and jS 1 : 

F Ax = Tcospcosa-D(V,a)&T-D(V,a) (9) 

F Az = -Tsma-L(V,a)^-L(V,a) (10) 

Fa t = -rsin/3cosa + 7 aero (F,/3)^7 aero (F,j3)(ll) 

with the following expansions for lift L , drag D and 
sideforce 7 aero : 

D(V, a) = qS (^Cj) 0 -\-Cj) a (X -\-Cj)^ 2 cc 2 ^j (12) 

L(V,a) = qS(C Lo +C La a ) (13) 

faero (V,P) = qS(c Yf) p ) (14) 

where the dynamic pressure q— 1/2 pV 2 . Currently, these 
conventional expansions are used in this approach. However, 
for future work, more elaborate expansions can be relied on, 
e.g. where drag D also depends on the absolute value of the 
sideslip angle |J3|. 

The corresponding numerical values are: m— 120- 10 2 kg, 
g = 9.81 m/s 2 , W = mg , p = 1.225 kg/m 2 (sea level), S = 
260 m 2 , C Lq = 1.0656, C La = 6.0723, C Dq = 0.1599, C Da = 
0.5035, C Da2 = 2.1175, C 7/3 =-1.6. 

In the perspective of reachability from stable and con- 
trollable trim conditions, the primary states of interest are 
airspeed V and flight path angle y. Considering time scale 
separation as presented in Fig. 2, the virtual inputs for the 
slow dynamics are roll angle <p, angle of attack a , sideslip 
angle j3 and thrust T. This framework and combining Eqs. 
(7)-(14) allows to define the general dynamics f in Eq. (1) 
by the following differential equation: 


V 

7 


-m v2 c Do -g sinr 

— f cos y 

cos a cos j8 

(cos (p sin a cos j8 — sin (p sin j8 ) y 


-t yl ( C Da a + C Dal a 2 
§£f(Ci 0 +Q, a a) cos <p _ 


0 

L— 5^/3 sin <p 

(15) 


Assuming small aerodynamic angles, as earlier, simplifies 
the differential equation: 


V 


'-^V 2 C Do -gsin r 

+ 

m 

.7 


— y cos y 


0 


7 + 


-S V 2 (cD a a+C Dal a 2 y 

- I^ V ( C Lo+ C L a a)coS(p _ 


0 

-^VCy^sinv 


(16) 


^ote that the allowed ranges in this specific example are set at a e 
[0°; 14.5°], j8 G [— 5°;5°], as defined later. The small angle assumption holds 
up to ±30°. 


This results in a Hamiltonian function with decoupled 
virtual inputs T, a and /3. Roll angle (p is not decoupled, but 
will be treated as a discretely gridded input. This decoupling 
significantly promotes computational efficiency, which is 
crucial for on-line applications. The Hamiltonian function 
becomes: 

//(p,x,u) = -T-pi^-V 2 C D a 2 + 

m 2m ° 

+ (p 2 C La cos (p-piVC Da ) a + 

P2 2m VCY P Sm ^ (17) 

where p are the co-states of the value function: p\ = ^ and 
P2= This Hamiltonian is linear in thrust T and sideslip 
angle /3 and quadratic in angle of attack a. This structure al- 
lows for an efficient optimization routine over the inputs. The 
trim envelope boundaries are: V m { n = 60m /s, F max — 100 m/s 
and y m i n / m ax = ±10°. The allowed ranges of the virtual inputs 
are: T min = 20546A, T max = 4109207V, a min = 0°, (W = 
14.5° (no stall), 0 min y max = ±60°, /3 min/max = ±5°. The 
maximizers T, a and /3 depend on the sign of the costates p\ 
and p 2 - Recall that V > 0 > 0,C^ 2 > 0,Ci a > 0 ,Cy^ < 
0. Due to the underlying physics, no sign changes for these 
parameters are to be expected in case of structural changes. 
Define p = P2CL(X ty*ycJ )l VCe>(x and a = amin + amax . Then the 
optimizing control inputs can be defined for invariance: 

• If p\ > 0 then t m T m [ n and 

if p > a then a = ( Xmm 

if p — a then a = ( Xmm or a max 

if p < a then a = ax 

• If Pi — 0 then T G [Tinin? ^max] and 

if p 2 > 0 then a = ( Xmm 
if P 2 = 0 then GC G [ofmm? <Anax] 
if P2 < 0 then a = a max 

• If p\ < 0 then f = T max and 

if P < tz m i n then 6c = (Xm i n 
if CCm\r\ ^ P ^ O^inax then CL — p 
if P ^ OCmax then CL — O-max 

• if p 2 sin (p > 0 then \ 3 = /3 m i n 

• if p 2 sin (p = 0 then j3 G [/3 min ;/3 max ] 

• if p 2 sin (p < 0 then \ 3 = /3 max 

For the purpose of maximizing the cost function with 
respect to the uncertainties A, the Hamiltonian from Eq. (17) 
can be rewritten, this time including parts independent of the 
inputs T, a or /3 but with some aerodynamic derivative(s): 

H{ p,x,u) = -p l ^V 2 C Do +p 2 ^VC Lo cos(p + 

~ P^V 2 {c Dl ,a + C Da2 a})+ (18) 

ps ps 

+ P 2 V-C La ^< ? a-p 2 -VCY !i sin <pfi 

It can be observed that the aerodynamic derivatives all appear 
linearly in an uncoupled way, which allows a similar proce- 
dure to solve the optimization as previously. By rewriting the 


Hamiltonian as a summation of terms, where each term is 
a multiplication of a variable involving a costate, a constant 
factor and a derivative, one can determine the sign of this 
factor, which consists of the predefined physical parameters: 

H (p,x,u) = -p l ^V 2 C Do +p 2 ^Vcos(pC Lo + 

>0 >0 
pS 9 9 pS 

~ P\^-V a C D 2 +p 2 —V a cos q>C La + 

2m ^ a 2m 

>o >o 

- P\ ^ v2ocCd « —pisintpP 2m VCr fi ( 19 ) 
>o x 

where it should be noted that <p G [— 60°;60°],a G 
[0°; 14.5°],J3 G [—5°; 5°]. Furthermore airspeed V > 0 and for 
the aerodynamic derivatives, it is known that Co a > 0 ,Cd 2 > 
0 ,Cl u > 0,Cyp < 0. Due to the underlying physics, no sign 
changes for these parameters are to be expected in case of 
uncertainty or structural changes. 

Based on this formulation, optimal control inputs for the 
aerodynamic derivatives can be defined as given in Table I 
where C. = C. max and C. = C. min . 

With this information, it is possible to create an entire 
“uncertainty band” around the envelope, however, here focus 
will be placed on the “worst-case” minimal size envelope. 

Fig. 4 compares the 3D envelopes with and without 
uncertainty, where two levels of uncertainty have been con- 
sidered here, namely 10% and 20% of uncertainty on all 
aerodynamic derivatives. For the purpose of this example, 
identical ratios of standard deviations over nominal values 
have been defined for all derivatives, but the algorithm is 
capable to deal with individual standard deviations which can 
vary between the different aerodynamic derivatives. It can be 
clearly seen in Fig. 4 that larger degrees of uncertainty result 
in more significant shrinking of the envelope, since this is a 
“worst-case” minimal size envelope. 

Fig. 5 analyzes the V, y maneuvering envelope for different 
values of bank angle <p, including robustness for uncertainties 
of 10% and 20%. By comparing Fig. 5(a), 5(b) and 5(c), 
it can be seen that larger bank angles have an influence 
on the climb capability of the aircraft. This is due to 
the physical principle that climb capability of lift force is 
provided through Lcos(p, which confirms a smaller decrease 
for smaller bank angles (up to (p = 25° as shown in Fig. 5(b)) 
but a much more significant change for larger bank angles 
as can be seen in Fig. 5(c). 

V. VALIDATION OF RESULTS 

The aforementioned results have been validated by means 
of Monte Carlo analyses. For this purpose, different bang- 
bang input signals have been inserted in the system. The 
extreme values of these signals correspond with the range 
limits. The time instant for the step change and the initial 
value for the input (maximum or minimum) vary over the 
Monte Carlo analysis. Running a nonlinear simulation of the 


TABLE I 

Optimal control inputs for robustness against uncertain aerodynamic derivatives 


sign of costate 


mimmizer 


maximizer 


Pi >0 
Pi< 0 

p 2 >0 

P2<0 

P 2 sin (p • /3 >0 


Cd 0 = Cd 0 , Co a = CDa , Cd u 2 = Co a 2 
Cd 0 = C Do , C Da = C Da , Co a2 = C Da2 
CLq = Cz 0 A a =Q La 
C L q =C Lq ,C L(X =C La 


C Y o = C Yr 


Cd 0 — C Dq , Co a — C Da , Co a 2 — C d ^ 2 

Cd {) = c Do , c D(y = c Da , c Da 2 = c D(x2 

Cl 0 =C Lo ,C La =C La 
Cl 0 = Ql 0 j Cl u = 

^ =&B 



P 2 sin 9 • j 8 <0 


CVfl = Cy 


Cyg = Cy 


(b) 9 = 25° 


V|mft| 

(c) (p = 60° 


Vflmfcl 

(a) <p = 0 ° 


Fig. 5. V,y maneuvering envelope of RCAM model for time horizon T = 2s for different bank angles and different uncertainty levels (0%, 10% and 20% 
uncertainty). Smaller envelope areas correspond to larger uncertainty bounds. 


Reachable region envelope of RCAM model at time = 0.00087s 



Fig. 4. Comparison of 3D envelopes with and without uncertainty: nominal 
(green), 10 % uncertainty (blue), 20 % uncertainty (yellow) 


aircraft model for the same time span as the time horizon 
T = 2s in the reachability analysis and plotting the traces 
in the envelopes, results in Fig. 6. For initial conditions 
xo within the backwards reachable set &{T = 2 s,K), it is 
always possible to find at least one admissible input u(-) 
which will bring part of the state trajectory 0(t,/,x,u,A) 
towards the end point at T = 2s inside the trim envelope K. 
On the other hand, from outside the backwards reachability 
set &{T = 2 s,K), it is impossible for the state trajectory 
0 (t,*,x,u, A) to reach the reference envelope K within T = 
2s, independent from which input u is applied. Many more 
Monte Carlo validations have been performed for different 
initial conditions xo, which all confirm the accuracy of the 
envelope in a similar way as shown here. Moreover, these 
Monte Carlo analyses have been based on the non-simplified 
aircraft model. As such, it has been demonstrated that the 
simplifying assumption made in Eq. (16) is acceptable and 
does not significantly perturb the results. 

VI. CONCLUSIONS 

In this paper, a computationally efficient algorithm for 
estimating the safe maneuvering envelope of damaged air- 
craft has been discussed. The algorithm performs a robust 
reachability analysis through an optimal control formulation 
while making use of time scale separation and taking into 
account uncertainties in the aerodynamic derivatives. The 
safe maneuvering envelope is defined as the cross section 




Fig. 6. Backwards reachability for RCAM model for time horizon T = 2s, 
including Monte Carlo Analysis 


between the forwards reachable and backwards reachable 
sets, which have been calculated starting from the stable 
trim envelope. Moreover, the backwards reachable set can 
be considered as the survivable maneuvering envelope, from 
where it is possible to bring the aircraft back to a safe trim 
condition after an upset due to damage, turbulence, a wake 
encounter etc. Results were found to be consistent with the 
underlying physical principles. This approach differs from 
others since it is physically inspired. This more transparent 
approach allows interpreting data in each step, and it is 
assumed that the physics based approach will therefore 
facilitate certification for future real life applications. 
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